TUTORIAL REGRESI POHON
Pada tutorial ini, akan dibahas tentang langkah-langkah untuk pemodelan pohon regresi dengan menggunakan software R. Tutorial ini berfokus pada penerapan pohon regresi untuk memprediksi harga rumah di King County, Washington, USA.
Sebelum memulai ini, diperlukan untuk menginstall beberapa package R yaitu:
library(dplyr)
library(viridis)
library(rpart)
library(mlr)
library(plotly)
library(leaflet)
library(rpart.plot)
library(parallelMap)
library(viridis)Import data dari file csv ke R
Pada tutorial ini, data yang akan digunakan adalah data tentang harga penjualan rumah di King County, WDataashington, USA. Data ini terdiri dari 18 kolom yaitu:
- price : Harga rumah
- bedrooms : Jumlah tempat tidur pada setiap rumah
- bathrooms : Jumlah kamar mandi pada setiap kamar
- waterfront : Rumah yang terdapat air di depannya
- view : Banyaknya Rumah pernah dilihat oleh calon pembeli
- condition : Kondisi Rumah (1-5), semakin besar semakin baik.
- grade : Peringkat Rumah berdasrkan sistem peringkat di King County
- yr_built : Tahun Rumah dibangun
- yr_renovated : Tahun Rumah direnovasi
- long : Garis Bujur
- lat : Garis lintang
- m2_living : Luas Rumah (meter persegi)
- m2_lot : Luas Tanah (meter persegi)
- floors : Total Lantai Rumah
- m2_above : Luas rumah tanpa basement (meter persegi)
- m2_basement : Luas basement (meter persegi)
- m2_living15 : Luas Rumah Setelah Renovasi (meter persegi)
- m2_lot15 : Luas Tanah Setelah Renovasi (meter persegi)
Langkah pertama yang harus dilakukan adalah menengimpor data dari file kc_house_data1.csv,
dta=read.csv(file = "D:/Job/E learning/Regression Tree/kc_house_data1.csv",sep = ",")
dta=dta%>%mutate(waterfront=as.factor(waterfront),grade=as.factor(grade),
view=as.factor(view),condition=as.factor(condition))
dta=dta%>%mutate(price=price/1000)
glimpse(dta)## Observations: 21,613
## Variables: 18
## $ price <dbl> 221.90, 538.00, 180.00, 604.00, 510.00, 1225.00, ...
## $ bedrooms <int> 3, 3, 2, 4, 3, 4, 3, 3, 3, 3, 3, 2, 3, 3, 5, 4, 3...
## $ bathrooms <dbl> 1.00, 2.25, 1.00, 3.00, 2.00, 4.50, 2.25, 1.50, 1...
## $ waterfront <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ view <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0...
## $ condition <fct> 3, 3, 3, 5, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3...
## $ grade <fct> 7, 7, 6, 7, 8, 11, 7, 7, 7, 7, 8, 7, 7, 7, 7, 9, ...
## $ yr_built <int> 1955, 1951, 1933, 1965, 1987, 2001, 1995, 1963, 1...
## $ yr_renovated <int> 0, 1991, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ long <dbl> -122.257, -122.319, -122.233, -122.393, -122.045,...
## $ lat <dbl> 47.5112, 47.7210, 47.7379, 47.5208, 47.6168, 47.6...
## $ m2_living <dbl> 109.62554, 238.76071, 71.53531, 182.08988, 156.07...
## $ m2_lot <dbl> 524.9020, 672.8035, 929.0300, 464.5150, 750.6562,...
## $ floors <dbl> 0.0929030, 0.1858060, 0.0929030, 0.0929030, 0.092...
## $ m2_above <dbl> 109.62554, 201.59951, 71.53531, 97.54815, 156.077...
## $ m2_basement <dbl> 0.00000, 37.16120, 0.00000, 84.54173, 0.00000, 14...
## $ m2_living15 <dbl> 124.49002, 157.00607, 252.69616, 126.34808, 167.2...
## $ m2_lot15 <dbl> 524.9020, 709.6860, 748.9840, 464.5150, 697.0512,...fillColor = "#FFA07A"
fillColor2 = "#F1C40F"Didalam tanda petik dua adalah alamat tempat file csv berada.
#Explorasi Data
Setelah data dimpor ke R, langkah selanjutnya adalah memahami data. Hal ini bisa dilakukan dengan melihat visualisasi data dalam bentuk grafik. Pertama-tama, akan dilihat grafik histogram dari price
plot_ly(data=dta,x=~price,type="histogram")Berdasarkan histogram, terdapat rumah-rumah yang sangat mahal pada data ini, seperti ada rumah yang seharga 7 juta dolar. Namun, sebagian besar dari rumah yang berada di King County, Washington, USA memiliki harga kurang dari 1 juta dollar. Setelah itu akan dilihat tentang bagaimana variabel-varibel selain price yang berperan sebagai variabel bebas memiliki keterkaitan dengan price.Variabel bebas pertama yang akan ditelusuri adalah variabel bedrooms.
dta%>%
group_by(bedrooms) %>%
summarise(PriceMean = mean(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(bedrooms = reorder(bedrooms,PriceMean)) %>%
arrange(desc(PriceMean)) %>%
ggplot(aes(x = bedrooms,y = round(PriceMean))) +
geom_bar(stat='identity',colour="white", fill = "#FFA07A") +
geom_text(aes(x = bedrooms, y = 1, label = paste0("(",round(PriceMean),")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'bedrooms',
y = 'Mean Price',
title = 'House Mean Price Based on Bedrooms') +
coord_flip() +
theme_bw()
Dari gambar diatas, dapat dilihat bahwa rata-rata harga rumah tidak dipengaruhi oleh banyaknya kamar. Sebagai contoh, rumah yang memliki kamar sebanyak 33 masih jauh lebih murah daripada rumah yang memiliki kamar sebanyak 5.
dta%>%
group_by(bathrooms) %>%
summarise(PriceMean = mean(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(bathrooms = reorder(bathrooms,PriceMean)) %>%
arrange(desc(PriceMean)) %>%
ggplot(aes(x = bathrooms,y = round(PriceMean))) +
geom_bar(stat='identity',colour="white", fill = "#FFA07A") +
geom_text(aes(x = bathrooms, y = 1, label = paste0("(",round(PriceMean),")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'bathrooms',
y = 'Mean Price',
title = 'House Mean Price Based on Bathrooms') +
coord_flip() +
theme_bw()
#Grade and Price
We examine how the Grade affects the price
dta %>%
group_by(grade) %>%
summarise(PriceMean = mean(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(grade = reorder(grade,PriceMean)) %>%
arrange(desc(PriceMean)) %>%
ggplot(aes(x = grade,y = round(PriceMean))) +
geom_bar(stat='identity',colour="white", fill = fillColor) +
geom_text(aes(x = grade, y = 1, label = paste0("(",round(PriceMean),")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'grade',
y = 'Mean Price',
title = 'grade and Mean Price') +
coord_flip() +
theme_bw()
waterfront and Price
We examine how the WaterFront affects the price
dta %>%
group_by(waterfront) %>%
summarise(PriceMedian = median(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(waterfront = reorder(waterfront,PriceMedian)) %>%
arrange(desc(PriceMedian)) %>%
ggplot(aes(x = waterfront,y = PriceMedian)) +
geom_bar(stat='identity',colour="white", fill = fillColor2) +
labs(x = 'waterfront',
y = 'Median Price',
title = 'waterfront and Median Price') +
theme_bw()
m2 Living and Price
We examine m2 Living and Price and plot a scatter plot
dta %>%
filter(!is.na(price)) %>%
filter(!is.na(m2_living)) %>%
ggplot(aes(x=m2_living,y=price))+
geom_point(color = "blue")+
stat_smooth(aes(x=m2_living,y=price),method="lm", color="red")+
theme_bw()+
theme(axis.title = element_text(size=16),axis.text = element_text(size=14))+
xlab("(m2 Living)")+
ylab("Price")
m2 Lot and Price
Plot 1
We examine m2 Lot and Price and plot a scatter plot
dta %>%
filter(!is.na(price)) %>%
filter(!is.na(m2_lot)) %>%
ggplot(aes(x=m2_lot,y=price))+
geom_point(color = "orange")+
stat_smooth(aes(x=m2_lot,y=price),method="lm", color="red")+
theme_bw()+
theme(axis.title = element_text(size=16),axis.text = element_text(size=14))+
xlab("(m2 Lot)")+
ylab("Price")
Plot 1 with Limits in X axis
We examine m2 Lot and Price and plot a scatter plot
dta %>%
filter(!is.na(price)) %>%
filter(!is.na(m2_lot)) %>%
ggplot(aes(x=m2_lot,y=price))+
geom_point(color = "orange")+
scale_x_continuous(limits=c(0,500e3)) +
stat_smooth(aes(x=m2_lot,y=price),method="lm", color="red")+
theme_bw()+
theme(axis.title = element_text(size=16),axis.text = element_text(size=14))+
xlab("(m2 Lot)")+
ylab("Price")
Year Built and Price
We examine the Year Built and Price and plot a bar plot.
dta %>%
group_by(yr_built) %>%
summarise(PriceMean = mean(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(yr_built = reorder(yr_built,PriceMean)) %>%
arrange(desc(PriceMean)) %>%
head(10) %>%
ggplot(aes(x = yr_built,y = round(PriceMean))) +
geom_bar(stat='identity',colour="white", fill = fillColor2) +
geom_text(aes(x = yr_built, y = 1, label = paste0("(",round(PriceMean),")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'year built',
y = 'Mean Price',
title = 'year built and Mean Price') +
coord_flip() +
theme_bw()
Year Renovated and Price
We examine the Year Renovated and Price and plot a bar plot.
dta %>%
group_by(yr_renovated) %>%
summarise(PriceMean = mean(price, na.rm = TRUE)) %>%
ungroup() %>%
mutate(yr_renovated = reorder(yr_renovated,PriceMean)) %>%
arrange(desc(PriceMean)) %>%
head(10) %>%
ggplot(aes(x = yr_renovated,y = round(PriceMean))) +
geom_bar(stat='identity',colour="white", fill = fillColor2) +
geom_text(aes(x = yr_renovated, y = 1, label = paste0("(",round(PriceMean),")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'year renovated',
y = 'Mean Price',
title = 'year renovated and Mean Price') +
coord_flip() +
theme_bw()
#Maps of Houses
Houses near the coast are costlier.
dta$PriceBin<-cut(dta$price, c(0,250,500,750,1000,2000,999000))
center_lon = median(dta$long,na.rm = TRUE)
center_lat = median(dta$lat,na.rm = TRUE)
factpal <- colorFactor(c("black","blue","yellow","orange","#0B5345","red"),
dta$PriceBin)
leaflet(dta) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal(PriceBin)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal, values = ~PriceBin,
title = "House Price Distribution",
opacity = 1)Price Bins Count
Most of the houses are in the range 250 thousand to 500 thousands. The next highest categories are
500 to 750 thousand
0 to 250 thousand
750 thousand to 1 million
1 million to 2 million
and the least is above 2 million
dta %>%
mutate(PriceBin = as.factor(PriceBin)) %>%
group_by(PriceBin) %>%
dplyr::summarise(Count = n()) %>%
ungroup() %>%
mutate(PriceBin = reorder(PriceBin,Count)) %>%
arrange(desc(Count)) %>%
ggplot(aes(x = PriceBin,y = Count)) +
geom_bar(stat='identity',colour="white", fill = fillColor2) +
geom_text(aes(x = PriceBin, y = 1, label = paste0("(",Count,")",sep="")),
hjust=0, vjust=.5, size = 4, colour = 'black',
fontface = 'bold') +
labs(x = 'PriceBin',
y = 'Count',
title = 'PriceBin and Count') +
coord_flip() +
theme_bw()
#Pemodelan pohon regresi Setelah kita memahami data, selanjutkan akan dilakukan pemodelan dengan menggunakan pohon regresi. Sebelum melakukan pemodelan, terlebih dahulu akan dilakukan pembagian data menjadi data training dan testing. Pembagian data ini dimaksudkan agar model yang dibagun tidak hanya dapat memprediksi data training saja, namun juga bisa memprediksi data diluar training dengan baik. Pada tutorial ini, pembagian data dilakukan dengan persentase 80% untuk data training dan 20% untuk data testing. Pada Software R, hal ini bisa dilakukan dengan sintaks dibawah ini dengan catatan harus menginstall terlebih dahulu menginstall package Caret
dta=dta%>%select(-PriceBin)
set.seed(123)
idx_train=caret::createDataPartition(dta$price,p=0.8,list=F)
train1=dta[idx_train,]
test=dta[-idx_train,]Salah satu hasil dari pohon regresi adalah grafik pohon. Grafik pohon digunakan untuk mengidentifikasi proses spliting yang terjadi antara variabel bebas pada pohon regresi.
mod=rpart(price~.,data=train1,method = "anova")
prp(mod,box.palette = viridis::viridis(n=4,alpha=0.5))
Berdasarkan grafik pohon diatas, ada beberapa hal yang bisa kita lihat. Pertama, pohon regresi untuk data ini dibangun dari 5 variabel bebas yaitu grade, lat, m2_living, yr_built dan long. Hal ini berarti 12 varibel bebas lainya dikeluarkan dari model pohon regresi. Root node pada gambar diatas adalah variabel grade, dimana jika peringkatnya 1,3,4,5,6,7,8 akan diteruskan ke sebelah kiri. Kata yes dan no menandakan kebenaran dari pernyataan sebelumnya. Jika bernilai benar maka pernyataan sebelumnya akan diteruskan ke kiri begitupun sebaliknya. Pernyataan yes dan no ini juga beraku pada node lainnya.Kemudian, untuk menjelaskan proses pohon regresi diatas akan digunakan salah satu alur dari atas ke bawah. Alur ini menjelaskan bagaimana domain (daerah) pada variabel bebas di split. Saat grade tidak bernilai 1,3,4,5,6,7,8, maka selanjutnya daerah dengan nilai grade 2,9,10,11,13 akan di split dengan variable m2_living. Kemudian saa nilai variabel m2_living kurang dari 377 maka akan di split lagi dengan variabel lat. Terakhir, saat daerah lat kurang dari 48 maka split akan berhenti. Ketika split berhenti, nilai prediksi akan dimunculkan. Secara singkat, bisa dikatakan bahwa model pohon regresi akan memprediksi price sebesar 552 saatgrade bernilai 2,9,10,11,13, m2_living kurang dari 377 dan lat kurang dari 48.
Setelah kita mengetahui bagaimana pohon regresi memodelkan data ini, selanjutnya akan dievaluasi bagaimana model ini dalam memprediksi data yang tidak dimasukan saat pemodelan (data test).
regtree_task=makeRegrTask(data = train1,target = "price")
regtree_learner=makeLearner('regr.rpart')
regtree_train=train(learner = regtree_learner,task = regtree_task)
regtree_pred=predict(regtree_train,newdata = test)
regtree_eval=performance(regtree_pred,measures = list(rmse,spearmanrho))Dari syntax diatas, akan menghasilkan output nilai RMSE dan spearmanrho yang digunakan untuk mengevaluasi hasil prediksi, berikut hasilnya
print(regtree_eval)## rmse spearmanrho
## 221.0285753 0.7857423Nilai RMSE (Root Mean Square Error) merupakan nilai rata-rata selisih antara price hasil prediksi dengan price sebenarnya. Sedangkan Spermanrho merupakan nilai korelasi yang dihitung menggunakan metode spearman. Nilai Spearmanrho ini berada pada selang 0 sampai 1. Semakin besar nilai spearmanrho maka hasil prediksi model pohon regresi semakin mendekati nilai price yang sebenarnya.RMSE sulit digunakan untuk menilai apakah pohon regresi memprediksi dengan baik atau tidak, karena nilai RMSE tidak memiliki batas minimum ataupun maksimum. RMSE akan lebih efektif digunakan jika terdapat lebih dari satu model. Sementara itu, Spermanrho bisa digunakan untuk menilai hasil prediksi dari model walaupun hanya terdapat satu model karena memiliki batas minimun dan maksimumnya.Berdasarkan hasil output diatas nilai spearmanrho adalah 0.785. Hal ini berarti prediksi pohon regresi sudah sangat baik. #Tuning hyperparameter Tuning hyperparameter dilakukan dengan tujuan meningkatkan akurasi model. hyperparameter adalah parameter yang terdapat pada model dimana nilainya ditentukan dengan trial and error. Hyperparamater pada pohon regresi adalah minimal split (minsplit), minimal bucket(minbucket), complexity parameter(cp), dan maximum depth(maxdepth). Minimal split merupakan parameter yang mengontrol banyaknya observasi minimum yang diperlukan pada suatu node sehingga split bisa dilakukan. Minimal bucket meruapakan parameter yang mengontrol banyaknya observasi minimum yang harus ada pada setiap node (biasanya bernilai sepertiga dari minsplit).Sebagai contoh, saat minsplit diatur bernilai 10 (minbucket akan bernilai 3) dengan data yang memiliki observasi 10 maka pohon regresi hanya akan memiliki root node. Hal ini dikarenakan walaupun berdasarkan minsplit root node bisa split namun karena node hasil split itu memiliki obeservasi 0 (dimana seharusnya observasi minimum pada setiap node sama dengan minbucket ) maka root node tidak bisa di split. Selanjutnya, complexity parameter merupakan parameter yang mengontrol spliting bisa dilakukan atau tidak. KOntrol dilakukan dengan apakah nilai MSE turun atau tidak dengan nilai minimum sebesar cp. Terakhir, maximum depth merupakan parameter untuk mengontorol kedalaman maksimum yang diperbolehkan pada pohon (pada plot bisa dilihat seberapa panjang gambar ke bawah), dengan root node dihitung sebagai depth 0.
Pada Software R, tuning hyperparameter bisa dilakukan dengan Grid Search, yaitu suatu metode pencarian hyperparameter yang optimal dengan mempertimbangkan nilai ukuran kebaikan model tertentu (seperti RMSE dan Spearmanrho).Langkah pertama yang dilakukan dalam Grid Search adalah user harus membuat list hyperparameter yang ingin dicobakan. Kemudian, secara otomatis grid search akan menampilkan hasil hyperparameter yang optimal.
regtree2_task=makeRegrTask(data = train1,target = "price")
regtree2_ps=makeParamSet(
makeDiscreteParam("minsplit",values=seq(10,50,length.out = 5)),
makeDiscreteParam("cp",values=seq(0.001,0.01,length.out = 5)),
makeDiscreteParam("maxdepth",values = seq(10,30,length.out = 5))
)Script R makeDiscreteParam adalah script untuk membuat list hyperparameter yang ingin dicobakan. Minimum bucket tidak dimasukan karena minbucket bergantung dari minsplit.
parallelStart(mode="socket",cpus = 4)## Starting parallelization in mode=socket with cpus=4.regtree2_ctrl=makeTuneControlGrid()
regtree2_rdesc=makeResampleDesc("CV",iter=10)
regtree2_tune=tuneParams("regr.rpart",
task = regtree2_task,par.set = regtree2_ps,
resampling =regtree2_rdesc,control = regtree2_ctrl,
measures = list(spearmanrho))## [Tune] Started tuning learner regr.rpart for parameter set:## Type len Def Constr Req Tunable
## minsplit discrete - - 10,20,30,40,50 - TRUE
## cp discrete - - 0.001,0.00325,0.0055,0.00775,0.01 - TRUE
## maxdepth discrete - - 10,15,20,25,30 - TRUE
## Trafo
## minsplit -
## cp -
## maxdepth -## With control class: TuneControlGrid## Imputation value: 1## Exporting objects to slaves for mode socket: .mlr.slave.options## Mapping in parallel: mode = socket; cpus = 4; elements = 125.## [Tune] Result: minsplit=20; cp=0.001; maxdepth=10 : spearmanrho.test.mean=0.8759661parallelStop()## Stopped parallelization. All cleaned up.Selanjutnya, script diatas merupakan script untuk menjalankan grid search dengan ukuran kebaikan model spearmanrho. Setelah nilai hyperparameter optimal ditemukan maka langkah selanjutnya adalah menerapkanya untuk memprediksi data testing,
regtree2_learner=makeLearner('regr.rpart',par.vals = regtree2_tune$x)
regtree2_train=train(learner = regtree2_learner,task = regtree2_task)
regtree2_pred=predict(regtree2_train,newdata = test)
regtree2_eval=performance(regtree2_pred,measures = list(rmse,spearmanrho))Kemudian, berikut adalah hasil rmse dan spearmanrho
print(regtree2_eval)## rmse spearmanrho
## 173.9054218 0.8733146Dapat dilihat bahwa spearmanrho dari pohon regresi yang hyperparameternya optimal mengalami peningkatan yang cukup besar, yaitu dari 0.785 menjadi 0.873. Kemudian, grafik pohon dari model ini bisa dilihat sebagai berikut:
mod1=rpart(price~.,data=train1,method = "anova",minsplit=20,cp=0.001,maxdepth=10)
prp(mod1,box.palette = viridis(n=4,alpha=0.5))
#Ringkasan 1. Regresi pohon memiliki empat hyperparameter, yaitu minsplit,minbucket,cp dan maxdepth 2. Tuning hyperparametric bisa meningkatan nilai akurasi model 3. Grafik pohon digunakan untuk mengidentifikasi proses spliting yang terjadi antara variabel bebas pada pohon regresi. 4. Korelasi Spearmanrho bisa digunakan untuk mengevaluasi kebaikan model walaupun hanya satu model saja